library(readxl)
library(httr)
library(arules)
library(tidyr)
library(dplyr)
library(lubridate)
library(stringr)
library(ggplot2)
library(purrr)
library(broom)
library(tibble)
library(plotly)
Odczyt danych przeprowdzony jest bezośrednio ze strony, aby zapewenić stałą wartość początkową dla każdego uruchomienia. Dodatkowo uzupełniane są brakujące dane identyfikatory pacjentów PATIENT_ID. Wartości brakujące wynikają ze sposobu przetwarzania scalonych komórek pliku źródłowego.
# https://stackoverflow.com/a/41368947
url <- "http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/wuhan_blood_sample_data_Jan_Feb_2020.xlsx"
GET(url, write_disk(tf <- tempfile(fileext = ".xlsx")))
df <- read_excel(tf) %>% fill(PATIENT_ID)
Surowy zbiór danych zawiera 81 atrybutów i 6120 wierszy.
Atrybuty zbioru danych można potwierdzić na dwie grupy:
Poza danymi niedostępnymi (NA) w zbiorze danych wsytępują dane brakujące, oznakowane wartością -1. Poniższa tabela pokazuje ilość takich wartości w każdej kolumnie, o ile te wartości występują.
df %>% pivot_longer(-c(RE_DATE,`Admission time`, `Discharge time`)) %>% filter(value==-1) %>% count(name)
## # A tibble: 2 x 2
## name n
## <chr> <int>
## 1 2019-nCoV nucleic acid detection 501
## 2 Platelet count 4
Dla kolumny Platelet count warkości brakujących jest mało, liczby -1 zostaną zamienione na NA.
df <- df %>% mutate(`Platelet count` = if(is.na(`Platelet count`) || `Platelet count` != -1) `Platelet count` else NA)
Dla kolumny 2019-nCoV nucleic acid detection sprawdzimy całkowitą ilość wartości różnych od NA.
df %>%
pivot_longer(-c(RE_DATE,`Admission time`, `Discharge time`)) %>%
filter(name=="2019-nCoV nucleic acid detection") %>%
filter(!is.na(value)) %>% count(name)
## # A tibble: 1 x 2
## name n
## <chr> <int>
## 1 2019-nCoV nucleic acid detection 501
Wszystkie wartości w kolumnie 2019-nCoV nucleic acid detection różne od NA wynoszą -1, w związku z tym ta kolumna jest właściwie nieprzydatna. Mimo to na ten moment postanowiono nie zmieniać wartości w tej kolumnie, ewentualne zmiany można wykonać w dalszych krokach.
Pierwszych siedem kolumn to atrybuty ogólne, reprezentują następująco:
Nazwy wyżej wymienionych atrybutów zostaną przeformatowane, będą pisane wielkimi literami, a ewentualne spacje będą zamieniane na znak podkreślinka _, w ten sposób atrbuty ogólne będą łatwo odróżnialne od pozostałych atrybutów.
df <- df %>% rename_with(~ toupper(gsub(" ", "_", .x, fixed = TRUE)), 1:7)
Pola GENDER i OUTCOME zostaną dodatkowo zmienione na typ factor - z typu numerycznego zamienione zostaną na zmienne nominalne. Do tego trzeba przeprowadzić identyfikację typu płci.
Wartości płci nie są jawnie opisane, w artykule źródłowym została natomiast podana proporcja pacjentów danej płci, co powinno umożliwić okreslenie właściwego przypisania.
knitr::kable(
df %>% group_by(PATIENT_ID) %>% summarise_all(last) %>% count(GENDER) %>% rename("count" = "n"),
align="lc"
)
| GENDER | count |
|---|---|
| 1 | 224 |
| 2 | 151 |
Wniosek:
12df <- df %>% mutate(across("GENDER", ~factor(., levels=c(1,2), labels=c("MALE", "FEMALE")))) %>%
mutate(across("OUTCOME", ~factor(., levels=c(0,1), labels=c("CURED","DECEASED"))))
Poniżej można zapoznać się z analizą wartości dla atrybutów ogólnych. Wartość atrybutu RE_DATE ma unikalną wartość dla kazdego wiersza.
knitr::kable(summary(select(df, RE_DATE)))
| RE_DATE | |
|---|---|
| Min. :2020-01-10 19:45:00 | |
| 1st Qu.:2020-02-04 13:44:00 | |
| Median :2020-02-09 12:42:30 | |
| Mean :2020-02-08 07:00:02 | |
| 3rd Qu.:2020-02-13 10:34:00 | |
| Max. :2020-02-18 17:49:00 | |
| NA’s :14 |
Pozostałe z atrybtów ogólnych, są wspólne dla każdego pacjenta, zastosowano transformację do zniwelowania skutków duplikacji.
options(knitr.kable.NA = '')
knitr::kable(
summary(
df %>%
select(c(1,3:7)) %>%
group_by(PATIENT_ID) %>%
summarise_all(last)
)
)
| PATIENT_ID | AGE | GENDER | ADMISSION_TIME | DISCHARGE_TIME | OUTCOME | |
|---|---|---|---|---|---|---|
| Min. : 1.0 | Min. :18.00 | MALE :224 | Min. :2020-01-10 15:52:20 | Min. :2020-01-23 09:09:23 | CURED :201 | |
| 1st Qu.: 94.5 | 1st Qu.:46.00 | FEMALE:151 | 1st Qu.:2020-02-01 19:27:40 | 1st Qu.:2020-02-11 13:39:21 | DECEASED:174 | |
| Median :188.0 | Median :62.00 | Median :2020-02-04 22:30:34 | Median :2020-02-16 17:40:07 | |||
| Mean :188.0 | Mean :58.83 | Mean :2020-02-04 20:13:51 | Mean :2020-02-15 16:42:59 | |||
| 3rd Qu.:281.5 | 3rd Qu.:70.00 | 3rd Qu.:2020-02-10 04:11:10 | 3rd Qu.:2020-02-19 11:47:14 | |||
| Max. :375.0 | Max. :95.00 | Max. :2020-02-17 21:30:07 | Max. :2020-03-04 16:21:51 |
Liczba unikatowych identyfikatorów pacjentów wynosi 375.
Ta sekcja poświęcona jest analizie długości hospitalizacji poszczególnych grup pacjentów.
| HOSPITALIZATION_TIME | quantile |
|---|---|
| 1 days | min |
| 6 days | 1st |
| 11 days | 2nd |
| 17 days | 3rd |
| 36 days | max |
| OUTCOME | GENDER | AGE | HOSPITALIZATION_TIME |
|---|---|---|---|
| CURED | MALE | 54.5 | 15.0 days |
| CURED | FEMALE | 45.0 | 15.0 days |
| DECEASED | MALE | 68.0 | 7.0 days |
| DECEASED | FEMALE | 72.0 | 6.5 days |
Pozostałe 74 kolumn zawiera wartości atrybutów, które reprezentują rodzaje miar zebranych podczas każdego badania.
Nazwy kolumn zostały poddane następującym transformacjom:
(%) zamieniono na suffix _percent(#) usunięto i myślnika - zamieniono na podkreślenie dolne _ ze względu na łatwość odwołania się do kolumnPoniższy blok kodu szczegółowo ukazuje zatosowane transfomrmacje, jednocześnie prezentowana są ostatenczne nazwy kolumn wyników badań.
df <- df %>%
rename("NT-proBNP" = "Amino-terminal brain natriuretic peptide precursor(NT-proBNP)") %>%
rename("2019-nCov-detection" = "2019-nCoV nucleic acid detection") %>%
rename("aPPT" = "Activation of partial thromboplastin time") %>%
rename("ALP" = "Alkaline phosphatase") %>%
rename("AAT" = "aspartate aminotransferase") %>%
rename("FDPs" = "Fibrin degradation products") %>%
rename("ALT" = "glutamic-pyruvic transaminase") %>%
rename("hs-CRP" = "High sensitivity C-reactive protein") %>%
rename("hs-cTnI" = "Hypersensitive cardiac troponinI") %>%
rename("HCV-abs-quant" = "HCV antibody quantification") %>%
rename("HIV-abs-quant" = "HIV antibody quantification") %>%
rename("INR" = "International standard ratio") %>%
rename("LDH" = "Lactate dehydrogenase") %>%
rename("MCH" = "mean corpuscular hemoglobin") %>%
rename("MCHC" = "mean corpuscular hemoglobin concentration") %>%
rename("MCV" = "mean corpuscular volume") %>%
rename("MPV" = "Mean platelet volume") %>%
rename("P-LCR" = "platelet large cell ratio") %>%
rename("PDW" = "PLT distribution width") %>%
rename("q-t-pallidum-abs" = "Quantification of Treponema pallidum antibodies") %>%
rename("RCDW-SD" = "RBC distribution width SD") %>%
rename("RBC_count" = "Red blood cell count") %>%
rename("RCDW" = "Red blood cell distribution width") %>%
rename("TNF-alfa" = "Tumor necrosis factorα") %>%
rename("WBC_count" = "White blood cell count") %>%
rename("gamma-GT" = "γ-glutamyl transpeptidase") %>%
rename_with(~str_c(str_replace(.x, "\\(%\\)", ""), "_percent"), contains("(%)")) %>%
rename_with(~str_replace(.x, "\\(#\\)", "")) %>%
rename_with(tolower, matches("^[[:upper:]]{1}[[:lower:]]+", ignore.case = FALSE)) %>%
rename_with(~str_replace_all(.x, "[ |-]", "_"))
str_sort(colnames(df[,-(1:7)]))
## [1] "2019_nCov_detection" "AAT" "albumin"
## [4] "ALP" "ALT" "antithrombin"
## [7] "aPPT" "basophil_count" "basophil_percent"
## [10] "calcium" "corrected_calcium" "creatinine"
## [13] "D_D_dimer" "direct_bilirubin" "eGFR"
## [16] "eosinophil_count" "eosinophils_percent" "ESR"
## [19] "FDPs" "ferritin" "fibrinogen"
## [22] "gamma_GT" "globulin" "glucose"
## [25] "HBsAg" "HCO3_" "HCV_abs_quant"
## [28] "hematocrit" "hemoglobin" "HIV_abs_quant"
## [31] "hs_CRP" "hs_cTnI" "indirect_bilirubin"
## [34] "INR" "interleukin_10" "interleukin_1β"
## [37] "interleukin_2_receptor" "interleukin_6" "interleukin_8"
## [40] "LDH" "lymphocyte_count" "lymphocyte_percent"
## [43] "MCH" "MCHC" "MCV"
## [46] "monocytes_count" "monocytes_percent" "MPV"
## [49] "neutrophils_count" "neutrophils_percent" "NT_proBNP"
## [52] "P_LCR" "PDW" "PH_value"
## [55] "platelet_count" "procalcitonin" "prothrombin_activity"
## [58] "prothrombin_time" "q_t_pallidum_abs" "RBC_count"
## [61] "RCDW" "RCDW_SD" "serum_chloride"
## [64] "serum_potassium" "serum_sodium" "thrombin_time"
## [67] "thrombocytocrit" "TNF_alfa" "total_bilirubin"
## [70] "total_cholesterol" "total_protein" "urea"
## [73] "uric_acid" "WBC_count"
Niniejsza sekcja jest poświęcona wstępnemu przetwarzniu zbiorów danych, ze szczególnym uwzględnieniem usuwania i scalania wybranych rekordów.
W powyższej głównej zauważono wystąpienia wartości brakujących dla kolumny RE_DATE. Szybkie sprawdzenie wykazuje, że dla wierszy, gdzie nie ma atrybut RE_DATE przyjmuje wartość ‘NA’, wartości wszystkich atrybutów pochodzących z badań krwi, również przyjmują wartość ‘NA’.
all(sapply(df[is.na(df$RE_DATE), -c(1,3:7)], is.na))
## [1] TRUE
Fakt, że wiersze te są właściwie puste, jest podstawą do ich usunięcia ze zbioru.
df <- df[!is.na(df$RE_DATE),]
Po tej zmianie w zbiorze pozostało 6106 wierszy i 361 unikalnych identyfikatorów pacjentów.
Wstępnę analiza danych sugeruje wskazuje, że wiele wartości atrybutów jest oznaczona wartością ‘NA’. Poniższy wykresy wskazuje, ilość wystąpień rekordów badań z danę liczbą wartości (różnych od NA).
Z powyższego histogramu można wyciągnąć kilka wniosków
W poprzedniej sekcji stwierdzono, że można przeanalizować współwystępowanie ze sobą poszczególnych atrybutów w rekordach. W związku z tym tabela wyników badań zostanie przekształcona do zbioru binarnych transakcji.
transactions <- sapply(df[,-(1:7)], function(x) !is.na(x))
Sprawdzona zostanie ilość unikalnych transakcji.
length(unique(apply(transactions, 1, function(x) which(x))))
## [1] 176
W zbiorze 6106 transakcji jest zaledwie 176 różnych typów.
Aby sprawdzić, czy atrybuty wsytępują w grupach, wyszukiwane są maksymalne domknięte zbiory częste, dla mimalnej wartości support równej 5%.
itemsets <- apriori(transactions, parameter = list(target="closed frequent itemsets", support=.05, minlen=1, maxlen=30, maxtime=60))
## Apriori
##
## Parameter specification:
## confidence minval smax arem aval originalSupport maxtime support minlen
## NA 0.1 1 none FALSE TRUE 60 0.05 1
## maxlen target ext
## 30 closed frequent itemsets TRUE
##
## Algorithmic control:
## filter tree heap memopt load sort verbose
## 0.1 TRUE TRUE FALSE TRUE 2 TRUE
##
## Absolute minimum support count: 305
##
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[74 item(s), 6106 transaction(s)] done [0.00s].
## sorting and recoding items ... [63 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 done [26.87s].
## filtering closed item sets ... done [15.64s].
## sorting transactions ... done [0.00s].
## writing ... [130 set(s)] done [0.34s].
## creating S4 object ... done [1.01s].
mc_itemsets <- itemsets[is.maximal(itemsets)]
inspect(sort(mc_itemsets, by = 'support'))
## items support transIdenticalToItemsets count
## [1] {hemoglobin,
## eosinophils_percent,
## basophil_percent,
## platelet_count,
## monocytes_percent,
## RCDW,
## neutrophils_percent,
## MCV,
## hematocrit,
## WBC_count,
## MCHC,
## lymphocyte_count,
## RBC_count,
## eosinophil_count,
## neutrophils_count,
## MPV,
## RCDW_SD,
## lymphocyte_percent,
## P_LCR,
## monocytes_count,
## PDW,
## basophil_count,
## MCH,
## thrombocytocrit} 0.13609564 0.1354406 831
## [2] {glucose} 0.12692434 0.0000000 775
## [3] {serum_chloride,
## ALP,
## albumin,
## total_bilirubin,
## indirect_bilirubin,
## total_protein,
## urea,
## corrected_calcium,
## serum_potassium,
## direct_bilirubin,
## total_cholesterol,
## AAT,
## uric_acid,
## HCO3_,
## calcium,
## LDH,
## globulin,
## gamma_GT,
## hs_CRP,
## serum_sodium,
## ALT,
## eGFR,
## creatinine} 0.11021946 0.1094006 673
## [4] {hs_cTnI} 0.08303308 0.0000000 507
## [5] {2019_nCov_detection} 0.08205044 0.0000000 501
## [6] {NT_proBNP} 0.07779234 0.0000000 475
## [7] {procalcitonin} 0.07517196 0.0000000 459
## [8] {PH_value} 0.06288896 0.0000000 384
## [9] {ESR} 0.06272519 0.0000000 383
## [10] {prothrombin_time,
## antithrombin,
## prothrombin_activity,
## fibrinogen,
## thrombin_time,
## D_D_dimer,
## FDPs,
## INR,
## aPPT} 0.05371765 0.0000000 328
Otrzymaliśmy 10 maksymalnych domkniętych zbiorów częstych, z czego 7 to zbiory jednolementowe. Pozostałe 3 zbiory (1, 3 i 10) zawierają wiele elementów. Na podstawie elementów składowych, zbiory wieolemelemnetowe zostały zakwalifikowane następująco:
Na podstawie uzyskanego wyniku można wnioskować, że ilość elementów w rekordzie może wynikać z procesu technologicznego analizowania próbek krwi. Sugeruje to, że możliwe będzie scalanie rekordów, np na podstawie dnia wykonania badania (jeśli jednemu pacjentowi przeprowadzono wiele badań w ciagu dnia).
TODO
class_labels <- as(items(sort(mc_itemsets, by="support")), "list")
class_labels <- sapply(class_labels, function(x) x[1])
class_labels[1] <- "Complete blood count"
class_labels[3] <- "Blood biochemistry"
class_labels[10] <- "Coagulation"
classes <- as(items(sort(mc_itemsets, by="support")), "matrix")
labeled <- df %>% rowwise() %>%
mutate(CLASS = paste(
which(
apply(
(matrix(
!is.na(c_across(-(1:7))),
nrow=nrow(classes),
ncol=ncol(classes),
byrow=TRUE,
dimnames=dimnames(classes)
) & classes) == classes,
1,
all)
), collapse=""),
.after=OUTCOME ) %>%
mutate(CLASS = if(CLASS != "" && as.integer(CLASS)>10) "MULTIPLE" else CLASS) %>%
mutate(CLASS = if(CLASS != "" && CLASS != "MULTIPLE") class_labels[as.integer(CLASS)] else CLASS) %>%
mutate(CLASS = str_replace(CLASS, "^$", "NONE")) %>%
mutate(CLASS = str_trim(CLASS)) %>%
mutate(CLASS = factor(CLASS)) %>%
rowwise() %>%
mutate(SIZE = sum(!is.na(c_across(-c(1:8)))), .after=CLASS) %>%
select(1:9) %>%
ungroup()
labeled %>% group_by(CLASS) %>% summarize(CNT = n(), .groups="drop") %>% arrange(desc(CNT))
## # A tibble: 12 x 2
## CLASS CNT
## <fct> <int>
## 1 NONE 1433
## 2 Complete blood count 816
## 3 glucose 584
## 4 Blood biochemistry 507
## 5 2019_nCov_detection 500
## 6 MULTIPLE 485
## 7 hs_cTnI 386
## 8 PH_value 375
## 9 ESR 368
## 10 Coagulation 300
## 11 NT_proBNP 190
## 12 procalcitonin 162
group.plot <- ggplot(
labeled,
aes(x=RE_DATE, y=CLASS, color=CLASS)
) +
geom_jitter(size=0.2) +
scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b %Y") +
theme(legend.position="none")
ggplotly(group.plot)
TODO Analiza ilości badań przeprowadzonych jednemu pacjentowi.
patient_last_results <- df %>%
group_by(PATIENT_ID) %>%
fill(everything()) %>%
summarise_all(last) %>%
rowwise() %>%
mutate(number_of_tests=sum(!is.na(c_across(-(1:7)))), .after=OUTCOME)
ggplot(
patient_last_results,
aes(number_of_tests)
) + geom_histogram(binwidth = 1)
Zauważamy, że część pacjentów ma niewiele przeprowadzonych rodzajów badań.
TODO analityczne wybranie tej wartości
* procent pełnych kolumn * procent pozostawionych pacjentów * wartość progu
veryfication_df <- function(n_tests, cols_cnt, th_seq) {
cls_vec <- c()
pts_vec <- c()
for(th in th_seq) {
tmp <- n_tests %>% filter(number_of_tests < cols_cnt*th)
dt <- n_tests %>% filter(!(PATIENT_ID %in% tmp$PATIENT_ID))
fl_cols <- mean(colSums(sapply(dt[,-(1:8)], is.na)) == 0)
pnts <- length(unique(dt$PATIENT_ID)) / length(unique(n_tests$PATIENT_ID))
cls_vec <- append(cls_vec, fl_cols)
pts_vec <- append(pts_vec, pnts)
}
data.frame(threshold=th_seq, full_cols=cls_vec, patients_left=pts_vec)
}
to_plot <- veryfication_df(patient_last_results, ncol(df)-8, seq(0.5, 1.0, by=.025))
ggplot(to_plot, aes(x=threshold)) +
geom_line(aes(y=full_cols, color="Full columns")) +
geom_line(aes(y=patients_left, color="Patients left")) +
labs(color="Legend") +
geom_vline(xintercept=0.675, linetype = "longdash")
TODO
W zbiorze danych pozostawimy tylko tych, którzy mają ponad 0.675 wykonanych rodzajów badań.
threshold <- 0.675
patients_to_be_removed <- patient_last_results %>%
filter(number_of_tests < (ncol(df)-7) * threshold) %>%
select(PATIENT_ID)
df <- df %>% filter(PATIENT_ID %!in% patients_to_be_removed$PATIENT_ID)
Po usunięciu rekordów pacjentów, u których przeprowadzono mało badań, pozostało 354 unikalnych pacjentów.
knitr::kable(
df %>%
group_by(PATIENT_ID) %>%
fill(everything()) %>%
summarise_all(last) %>%
ungroup() %>%
select(-(1:7)) %>%
pivot_longer(everything()) %>%
mutate(value=is.na(value)) %>%
group_by(name) %>%
summarize(value=sum(value), .groups="drop") %>%
arrange(value)
)
| name | value |
|---|---|
| AAT | 0 |
| albumin | 0 |
| ALP | 0 |
| ALT | 0 |
| basophil_count | 0 |
| basophil_percent | 0 |
| creatinine | 0 |
| direct_bilirubin | 0 |
| eGFR | 0 |
| eosinophil_count | 0 |
| eosinophils_percent | 0 |
| gamma_GT | 0 |
| globulin | 0 |
| HCO3_ | 0 |
| hematocrit | 0 |
| hemoglobin | 0 |
| LDH | 0 |
| lymphocyte_count | 0 |
| lymphocyte_percent | 0 |
| MCH | 0 |
| MCHC | 0 |
| MCV | 0 |
| monocytes_count | 0 |
| monocytes_percent | 0 |
| neutrophils_count | 0 |
| neutrophils_percent | 0 |
| platelet_count | 0 |
| RBC_count | 0 |
| total_bilirubin | 0 |
| total_cholesterol | 0 |
| total_protein | 0 |
| urea | 0 |
| uric_acid | 0 |
| WBC_count | 0 |
| indirect_bilirubin | 1 |
| calcium | 2 |
| serum_chloride | 2 |
| serum_potassium | 2 |
| serum_sodium | 2 |
| corrected_calcium | 3 |
| hs_CRP | 3 |
| glucose | 5 |
| INR | 5 |
| prothrombin_activity | 5 |
| prothrombin_time | 5 |
| RCDW | 6 |
| RCDW_SD | 6 |
| MPV | 10 |
| P_LCR | 10 |
| PDW | 10 |
| thrombocytocrit | 10 |
| D_D_dimer | 14 |
| procalcitonin | 44 |
| aPPT | 59 |
| fibrinogen | 59 |
| thrombin_time | 59 |
| ESR | 69 |
| hs_cTnI | 69 |
| HBsAg | 82 |
| HCV_abs_quant | 82 |
| q_t_pallidum_abs | 82 |
| HIV_abs_quant | 83 |
| NT_proBNP | 90 |
| PH_value | 126 |
| interleukin_6 | 138 |
| 2019_nCov_detection | 140 |
| interleukin_1β | 140 |
| interleukin_2_receptor | 140 |
| interleukin_8 | 140 |
| TNF_alfa | 140 |
| interleukin_10 | 141 |
| ferritin | 143 |
| antithrombin | 153 |
| FDPs | 153 |
df <- df %>%
mutate(DAYS_TO_OUTCOME = as.numeric(difftime(
floor_date(DISCHARGE_TIME, "1 day"),
floor_date(RE_DATE, "1 day"), unit="days"
), unit="days"), .before=OUTCOME)
df %>% summarise(DAYS_TO_OUTCOME = quantile(DAYS_TO_OUTCOME, c(0, 0.25, 0.5, 0.75,1)), quantile = c("min", "1st", "2nd", "3rd", "max"))
## # A tibble: 5 x 2
## DAYS_TO_OUTCOME quantile
## <dbl> <chr>
## 1 -2 min
## 2 3 1st
## 3 7 2nd
## 4 13 3rd
## 5 35 max
ggplot(df, aes(DAYS_TO_OUTCOME)) + geom_histogram(binwidth = 1)
df <- df %>%
mutate(RE_DATE = floor_date(RE_DATE, '1 day')) %>%
group_by(PATIENT_ID, RE_DATE) %>%
fill(everything()) %>%
summarise_all(last) %>%
ungroup() %>%
rowwise() %>%
mutate(number_of_tests=sum(!is.na(c_across(-(1:8)))), .after=OUTCOME)
ggplot(df, aes(number_of_tests)) + geom_histogram(binwidth = 1)
df <- df %>% select(-number_of_tests)
ggplot(
df %>% pivot_longer(hs_cTnI:creatinine) %>% filter(!is.na(value)),
aes(x=OUTCOME, y=value)
) + geom_boxplot() +
scale_x_discrete(labels = abbreviate) +
facet_wrap(name ~ ., scales="free",ncol=5)
TODO
cor_mat <- cor(
x = df %>%
mutate(isMALE = as.numeric(GENDER=="MALE"),
isCURED = as.numeric(OUTCOME=="CURED"),
has_nCOV_test = as.numeric(!is.na(`2019_nCov_detection`))) %>%
select(-c(1:8), -"2019_nCov_detection"),
use="pairwise.complete.obs"
)
cor_df = data.frame(round(cor_mat,2)) %>%
rownames_to_column() %>%
pivot_longer(-rowname, names_to="colname")
cor_plot <- ggplot(cor_df, aes(colname, rowname, fill=value)) +
geom_tile() +
scale_fill_gradient2() +
theme(axis.text.x=element_text(angle = 90, hjust = 0))
ggplotly(cor_plot)
knitr::kable(
cor_df %>% filter(colname>rowname) %>% arrange(desc(abs(value))) %>% head(30)
)
| rowname | colname | value |
|---|---|---|
| MPV | P_LCR | 0.99 |
| INR | prothrombin_time | 0.99 |
| platelet_count | thrombocytocrit | 0.98 |
| HBsAg | interleukin_6 | 0.98 |
| direct_bilirubin | total_bilirubin | 0.98 |
| lymphocyte_percent | neutrophils_percent | -0.98 |
| procalcitonin | q_t_pallidum_abs | 0.95 |
| D_D_dimer | FDPs | 0.95 |
| P_LCR | PDW | 0.95 |
| MPV | PDW | 0.94 |
| ALT | interleukin_1β | 0.93 |
| lymphocyte_count | monocytes_count | 0.88 |
| serum_chloride | serum_sodium | 0.86 |
| eosinophil_count | eosinophils_percent | 0.86 |
| AAT | interleukin_1β | 0.86 |
| hematocrit | hemoglobin | 0.84 |
| MCH | MCV | 0.82 |
| indirect_bilirubin | total_bilirubin | 0.80 |
| interleukin_6 | interleukin_8 | 0.80 |
| RCDW | RCDW_SD | 0.79 |
| aPPT | prothrombin_time | 0.79 |
| monocytes_percent | neutrophils_percent | -0.77 |
| eGFR | urea | -0.77 |
| creatinine | urea | 0.75 |
| albumin | calcium | 0.74 |
| isCURED | neutrophils_percent | -0.73 |
| isCURED | lymphocyte_percent | 0.72 |
| neutrophils_count | neutrophils_percent | 0.71 |
| albumin | total_protein | 0.70 |
| lymphocyte_percent | neutrophils_count | -0.70 |
knitr::kable(
tidy(cor.test(df$LDH, as.numeric(df$OUTCOME=="CURED")))
)
| estimate | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|
| -0.6270327 | -24.33491 | 0 | 914 | -0.6648056 | -0.5860615 | Pearson’s product-moment correlation | two.sided |
correlations <- df %>% mutate(isCURED=as.numeric(OUTCOME=="CURED")) %>% select(-c(1:2,4:8)) %>%
pivot_longer(!isCURED, names_to="attribute", values_to="x_val") %>%
filter(!is.na(x_val)) %>%
nest(data = c(x_val, isCURED)) %>%
mutate(cor_test = map(data, ~tidy(cor.test(.x$x_val, .x$isCURED)))) %>%
unnest(cor_test) %>%
select(attribute, estimate, p.value, conf.low, conf.high) %>%
arrange(desc(abs(estimate))) %>% head(15)
## Warning: Problem with `mutate()` input `cor_test`.
## ℹ the standard deviation is zero
## ℹ Input `cor_test` is `map(data, ~tidy(cor.test(.x$x_val, .x$isCURED)))`.
## Warning in cor(x, y): the standard deviation is zero
knitr::kable(
correlations
)
| attribute | estimate | p.value | conf.low | conf.high |
|---|---|---|---|---|
| neutrophils_percent | -0.7301461 | 0 | -0.7586392 | -0.6988654 |
| lymphocyte_percent | 0.7234883 | 0 | 0.6915912 | 0.7525691 |
| albumin | 0.6880859 | 0 | 0.6524294 | 0.7207028 |
| prothrombin_activity | 0.6547754 | 0 | 0.6084767 | 0.6966318 |
| hs_CRP | -0.6509896 | 0 | -0.6910468 | -0.6069460 |
| D_D_dimer | -0.6376192 | 0 | -0.6819778 | -0.5885867 |
| LDH | -0.6270327 | 0 | -0.6648056 | -0.5860615 |
| neutrophils_count | -0.6083675 | 0 | -0.6470960 | -0.5665074 |
| FDPs | -0.6042409 | 0 | -0.6689583 | -0.5304310 |
| AGE | -0.5761382 | 0 | -0.6071595 | -0.5433633 |
| calcium | 0.5703332 | 0 | 0.5257539 | 0.6117881 |
| platelet_count | 0.5419441 | 0 | 0.4952125 | 0.5855487 |
| monocytes_percent | 0.5415337 | 0 | 0.4947996 | 0.5851444 |
| eGFR | 0.4909936 | 0 | 0.4402769 | 0.5385870 |
| urea | -0.4841428 | 0 | -0.5321758 | -0.4330031 |
tseries <- df %>%
select(1:8) %>%
group_by(PATIENT_ID) %>%
summarize(across(everything(), last), .groups="drop") %>%
select(-RE_DATE, -AGE, -GENDER, -DAYS_TO_OUTCOME, -PATIENT_ID) %>%
pivot_longer(ADMISSION_TIME:DISCHARGE_TIME, names_to="TYPE", values_to="DATE") %>%
mutate(DATE = floor_date(DATE, "days")) %>%
mutate(TYPE = str_replace(TYPE, "_TIME", "")) %>%
count(TYPE, DATE, OUTCOME) %>%
rename("COUNT" = "n")
tplot <- ggplot(tseries, aes(x=DATE, y=COUNT, fill=OUTCOME)) +
geom_col() +
scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b %Y") +
facet_grid(TYPE ~ .)
ggplotly(tplot)
unique(c(.packages(), loadedNamespaces()))
## [1] "plotly" "tibble" "broom" "purrr" "ggplot2"
## [6] "stringr" "lubridate" "dplyr" "tidyr" "arules"
## [11] "Matrix" "httr" "readxl" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
## [21] "Rcpp" "highr" "cellranger" "compiler" "pillar"
## [26] "tools" "digest" "viridisLite" "jsonlite" "evaluate"
## [31] "lifecycle" "gtable" "lattice" "pkgconfig" "rlang"
## [36] "cli" "crosstalk" "curl" "yaml" "xfun"
## [41] "withr" "knitr" "htmlwidgets" "generics" "vctrs"
## [46] "grid" "tidyselect" "data.table" "glue" "R6"
## [51] "fansi" "rmarkdown" "farver" "magrittr" "backports"
## [56] "scales" "htmltools" "ellipsis" "assertthat" "colorspace"
## [61] "labeling" "utf8" "stringi" "lazyeval" "munsell"
## [66] "crayon"